home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1996 September / macformat-041.iso / mac / Shareware City / Graphics / MacSPD / Sources / libtx.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-04-04  |  12.4 KB  |  520 lines  |  [TEXT/MMCC]
  1. /*
  2.  * libtx.c - a library of transformation tracking routines.
  3.  *
  4.  * Author:  Alexander Enzmann
  5.  *
  6.  */
  7.  
  8.  
  9. /*-----------------------------------------------------------------*/
  10. /* include section */
  11.  
  12. #include <stdio.h>
  13. #include <stdlib.h>
  14. #include <math.h>
  15. #include <string.h>
  16. #include "lib.h"
  17. #include "drv.h"
  18.  
  19.  
  20. /*-----------------------------------------------------------------*/
  21. /* defines/constants section */
  22.  
  23. typedef struct tx_struct *tx_ptr;
  24. typedef struct tx_struct {
  25.    MATRIX tx;
  26.    tx_ptr next;
  27.    };
  28.  
  29.  
  30. /*-----------------------------------------------------------------*/
  31. MATRIX IdentityTx =
  32.    {{1, 0, 0, 0},
  33.     {0, 1, 0, 0},
  34.     {0, 0, 1, 0},
  35.     {0, 0, 0, 1}};
  36. static MATRIX CurrentTx =
  37.    {{1, 0, 0, 0},
  38.     {0, 1, 0, 0},
  39.     {0, 0, 1, 0},
  40.     {0, 0, 0, 1}};
  41. static tx_ptr TxStack = NULL;
  42.  
  43. /* Return 1 if there is an active transformation, 0 if not */
  44. int
  45. lib_tx_active()
  46. {
  47.    int i, j;
  48.    for (i=0;i<4;i++)
  49.       for (j=0;j<4;j++)
  50.      if (fabs(CurrentTx[i][j] - (i == j ? 1.0 : 0.0)) > EPSILON)
  51.         return 1;
  52.    return 0;
  53. }
  54.  
  55. /* Copy the current transform into mat */
  56. void
  57. lib_get_current_tx(mat)
  58.    MATRIX mat;
  59. {
  60.    memcpy(mat, CurrentTx, sizeof(MATRIX));
  61. }
  62.  
  63. /* Copy matrix mat into the current transform */
  64. void
  65. lib_set_current_tx(mat)
  66.    MATRIX mat;
  67. {
  68.    memcpy(CurrentTx, mat, sizeof(MATRIX));
  69. }
  70.  
  71. #ifdef _DEBUG
  72. static void
  73. show_matrix(mat)
  74.    MATRIX mat;
  75. {
  76.    int i, j;
  77.    for (i=0;i<4;i++) {
  78.       for (j=0;j<4;j++)
  79.      fprintf(stderr, "%7.4g ", mat[i][j]);
  80.       fprintf(stderr, "\n");
  81.       }
  82. }
  83.  
  84. /* Show what the transform looks like */
  85. static void
  86. lib_tx_status()
  87. {
  88.    double trans[16];
  89.  
  90. #if 0
  91. fprintf(stderr, "Current tx:\n");
  92. show_matrix(CurrentTx);
  93. #endif
  94.  
  95.    if (lib_tx_active()) {
  96.       if (lib_tx_unwind(CurrentTx, trans)) {
  97.      fprintf(stderr, "Tx: scale     <%g,%g,%g>\n",
  98.          trans[U_SCALEX], trans[U_SCALEY], trans[U_SCALEZ]);
  99.      fprintf(stderr, "    rotate    <%g,%g,%g>\n",
  100.          trans[U_ROTATEX], trans[U_ROTATEY], trans[U_ROTATEZ]);
  101.      fprintf(stderr, "    translate <%g,%g,%g>\n",
  102.          trans[U_TRANSX], trans[U_TRANSY], trans[U_TRANSZ]);
  103.      }
  104.       else
  105.      fprintf(stderr, "Tx: failed to unwind\n");
  106.       }
  107.    else
  108.       fprintf(stderr, "Tx: Identity\n");
  109. }
  110. #endif
  111.  
  112. /* Turn a transformation matrix into the appropriate sequence of
  113.    translate/rotate/scale instructions */
  114. void
  115. lib_output_tx_sequence()
  116. {
  117.    MATRIX txmat;
  118.    double trans[16];
  119.    int i, j, tflag, rflag, sflag;
  120.  
  121.     lib_get_current_tx(txmat);
  122.  
  123.     if (!lib_tx_unwind(txmat, trans)) {
  124.      fprintf(stderr, "Tx: failed to unwind\n");
  125.      return;
  126.      }
  127.     if (fabs(trans[U_SCALEX] - 1.0) > EPSILON ||
  128.     fabs(trans[U_SCALEY] - 1.0) > EPSILON ||
  129.     fabs(trans[U_SCALEZ] - 1.0) > EPSILON)
  130.     sflag = 1;
  131.     else
  132.     sflag = 0;
  133.     if (fabs(trans[U_ROTATEX]) > EPSILON ||
  134.     fabs(trans[U_ROTATEY]) > EPSILON ||
  135.     fabs(trans[U_ROTATEZ]) > EPSILON) {
  136.     trans[U_ROTATEX] = RAD2DEG(trans[U_ROTATEX]);
  137.     trans[U_ROTATEY] = RAD2DEG(trans[U_ROTATEY]);
  138.     trans[U_ROTATEZ] = RAD2DEG(trans[U_ROTATEZ]);
  139.     rflag = 1;
  140.     }
  141.     else
  142.     rflag = 0;
  143.     if (fabs(trans[U_TRANSX]) > EPSILON ||
  144.     fabs(trans[U_TRANSY]) > EPSILON ||
  145.     fabs(trans[U_TRANSZ]) > EPSILON)
  146.     tflag = 1;
  147.     else
  148.     tflag = 0;
  149.     if (!sflag && !rflag && !tflag)
  150.        return;
  151.     switch (gRT_out_format) {
  152.     case OUTPUT_VIDEO:
  153.     case OUTPUT_DELAYED:
  154.     case OUTPUT_DXF:
  155.     case OUTPUT_PLG:
  156.     case OUTPUT_NFF:
  157.     case OUTPUT_OBJ:
  158.     case OUTPUT_QRT:
  159.     case OUTPUT_RAWTRI:
  160.         /* Can't do inline transforms in these renderers, the
  161.            code does the transformations on the shapes
  162.            themselves. */
  163.         break;
  164.  
  165.     case OUTPUT_RTRACE:
  166.         fprintf(gOutfile, "65 %ld ", gObject_count+1);
  167.         for (i=0;i<4;i++)
  168.                 for (j=0;j<4;j++)
  169.                     fprintf(gOutfile, "%g ", txmat[j][i]);
  170.             fprintf(gOutfile, "\n");
  171.     break;
  172.  
  173.     case OUTPUT_RWX:
  174.         /* Transforms are possible in RWX,
  175.            this code needs to be finished... */
  176.         if (tflag) {
  177.             tab_indent();
  178.             fprintf(gOutfile, "Translate %g %g %g\n",
  179.                 trans[U_TRANSX], trans[U_TRANSY], trans[U_TRANSZ]);
  180.             }
  181.         if (rflag) {
  182.             tab_indent();
  183.             fprintf(gOutfile, "Rotate 0 0 1 %g\n",
  184.                 trans[U_ROTATEZ]);
  185.             tab_indent();
  186.             fprintf(gOutfile, "Rotate 0 1 0 %g\n",
  187.                 trans[U_ROTATEY]);
  188.             tab_indent();
  189.             fprintf(gOutfile, "Rotate 1 0 0 %g\n",
  190.                 trans[U_ROTATEX]);
  191.             }
  192.         if (sflag) {
  193.             tab_indent();
  194.             fprintf(gOutfile, "Scale %g %g %g\n",
  195.                 trans[U_SCALEX], trans[U_SCALEY], trans[U_SCALEZ]);
  196.             }
  197.         break;
  198.  
  199.     case OUTPUT_RIB:
  200.         if (sflag) {
  201.             tab_indent();
  202.             fprintf(gOutfile, "Scale %#g %#g %#g\n",
  203.                 trans[U_SCALEX], trans[U_SCALEY], trans[U_SCALEZ]);
  204.             }
  205.         if (rflag) {
  206.             tab_indent();
  207.             fprintf(gOutfile, "Rotate %#g 1 0 0\n",
  208.                 trans[U_ROTATEX]);
  209.             tab_indent();
  210.             fprintf(gOutfile, "Rotate %#g 0 1 0\n",
  211.                 trans[U_ROTATEY]);
  212.             tab_indent();
  213.             fprintf(gOutfile, "Rotate %#g 0 0 1\n",
  214.                 trans[U_ROTATEZ]);
  215.             }
  216.         if (tflag) {
  217.             tab_indent();
  218.             fprintf(gOutfile, "Translate %#g %#g %#g\n",
  219.                 trans[U_TRANSX], trans[U_TRANSY], trans[U_TRANSZ]);
  220.             }
  221.         break;
  222.  
  223.     case OUTPUT_VIVID:
  224.         tab_inc();
  225.         if (sflag) {
  226.             tab_indent();
  227.             fprintf(gOutfile, "scale %g\n", trans[U_SCALEX]);
  228.             }
  229.         tab_dec();
  230.         if (rflag) {
  231.             tab_indent();
  232.             fprintf(gOutfile, "rotate %g %g %g\n",
  233.                 trans[U_ROTATEX], trans[U_ROTATEY], trans[U_ROTATEZ]);
  234.             }
  235.         if (tflag) {
  236.             tab_indent();
  237.             fprintf(gOutfile, "translate %g %g %g\n",
  238.                 trans[U_TRANSX], trans[U_TRANSY], trans[U_TRANSZ]);
  239.             }
  240.         break;
  241.  
  242.     case OUTPUT_RAYSHADE:
  243.         if (sflag) {
  244.             fprintf(gOutfile, " scale %g %g %g",
  245.                 trans[U_SCALEX], trans[U_SCALEY], trans[U_SCALEZ]);
  246.         }
  247.         if (rflag) {
  248.             fprintf(gOutfile, " rotate 1 0 0 %g",
  249.                 trans[U_ROTATEX]);
  250.             tab_indent();
  251.             fprintf(gOutfile, " rotate 0 1 0 %g",
  252.                 trans[U_ROTATEY]);
  253.             tab_indent();
  254.             fprintf(gOutfile, " rotate 0 0 1 %g",
  255.                 trans[U_ROTATEZ]);
  256.         }
  257.         if (tflag) {
  258.             fprintf(gOutfile, " translate %g %g %g",
  259.                 trans[U_TRANSX], trans[U_TRANSY], trans[U_TRANSZ]);
  260.         }
  261.         break;
  262.  
  263.     case OUTPUT_POVRAY_10:
  264.         if (sflag) {
  265.             fprintf(gOutfile, " scale %g %g %g",
  266.                 trans[U_SCALEX], trans[U_SCALEY], trans[U_SCALEZ]);
  267.         }
  268.         if (rflag) {
  269.             fprintf(gOutfile, " rotate %g %g %g",
  270.                 trans[U_ROTATEX], trans[U_ROTATEY], trans[U_ROTATEZ]);
  271.         }
  272.         if (tflag) {
  273.             fprintf(gOutfile, " translate %g %g %g",
  274.                 trans[U_TRANSX], trans[U_TRANSY], trans[U_TRANSZ]);
  275.         }
  276.         break;
  277.  
  278.     case OUTPUT_POVRAY_20:
  279.     case OUTPUT_POLYRAY:
  280.         if (sflag) {
  281.             fprintf(gOutfile, " scale <%g, %g, %g>",
  282.                 trans[U_SCALEX], trans[U_SCALEY], trans[U_SCALEZ]);
  283.         }
  284.         if (rflag) {
  285.             tab_indent();
  286.             fprintf(gOutfile, " rotate <%g, %g, %g>",
  287.                 trans[U_ROTATEX], trans[U_ROTATEY], trans[U_ROTATEZ]);
  288.         }
  289.         if (tflag) {
  290.             tab_indent();
  291.             fprintf(gOutfile, " translate <%g, %g, %g>",
  292.                 trans[U_TRANSX], trans[U_TRANSY], trans[U_TRANSZ]);
  293.         }
  294.         break;
  295.  
  296.     case OUTPUT_ART:
  297.         if (sflag) {
  298.             tab_indent();
  299.             fprintf(gOutfile, "scale(%g, %g, %g)\n",
  300.                 trans[U_SCALEX], trans[U_SCALEY], trans[U_SCALEZ]);
  301.         }
  302.         if (rflag) {
  303.             tab_indent();
  304.             fprintf(gOutfile, "rotate(%g, x)\n",
  305.                 trans[U_ROTATEX]);
  306.             tab_indent();
  307.             fprintf(gOutfile, "rotate(%g, y)\n",
  308.                 trans[U_ROTATEY]);
  309.             tab_indent();
  310.             fprintf(gOutfile, "rotate(%g, z)\n",
  311.                 trans[U_ROTATEZ]);
  312.         }
  313.         if (tflag) {
  314.             tab_indent();
  315.             fprintf(gOutfile, "translate(%g, %g, %g)\n",
  316.                 trans[U_TRANSX], trans[U_TRANSY], trans[U_TRANSZ]);
  317.         }
  318.         break;
  319.     }
  320. } /* lib_output_tx_sequence */
  321.  
  322. void
  323. lib_tx_pop()
  324. {
  325.    tx_ptr last_tx;
  326.  
  327.    if (TxStack == NULL) {
  328.       fprintf(stderr, "Attempt to pop beyond bottom of transform stack\n");
  329.       }
  330.    else {
  331.       last_tx = TxStack;
  332.       lib_copy_matrix(CurrentTx, last_tx->tx);
  333.       TxStack = TxStack->next;
  334.       free(last_tx);
  335.       }
  336. }
  337.  
  338. void
  339. lib_tx_push()
  340. {
  341.    tx_ptr new_tx;
  342.  
  343.     if ((new_tx = malloc(sizeof(struct tx_struct))) == NULL) {
  344.     fprintf(stderr, "Failed to allocate polygon data\n");
  345.     exit(EXIT_FAIL);
  346.     }
  347.     lib_copy_matrix(new_tx->tx, CurrentTx);
  348.     new_tx->next = TxStack;
  349.     TxStack = new_tx;
  350. }
  351.  
  352. void
  353. lib_tx_rotate(axis, angle)
  354.     int axis;
  355.     double angle;
  356. {
  357.     MATRIX mx1, mx2;
  358.  
  359.     lib_create_rotate_matrix(mx1, axis, angle);
  360.     lib_copy_matrix(mx2, CurrentTx);
  361.     lib_matrix_multiply(CurrentTx, mx1, mx2);
  362. }
  363.  
  364. void
  365. lib_tx_scale(vec)
  366.     COORD3 vec;
  367. {
  368.     MATRIX mx1, mx2;
  369.  
  370.     lib_create_scale_matrix(mx1, vec);
  371.     lib_copy_matrix(mx2, CurrentTx);
  372.     lib_matrix_multiply(CurrentTx, mx1, mx2);
  373. }
  374.  
  375. void
  376. lib_tx_translate(vec)
  377.     COORD3 vec;
  378. {
  379.     MATRIX mx1, mx2;
  380.  
  381.     lib_create_translate_matrix(mx1, vec);
  382.     lib_copy_matrix(mx2, CurrentTx);
  383.     lib_matrix_multiply(CurrentTx, mx1, mx2);
  384. }
  385.  
  386. /* From Graphics Gems II, "unmatrix" written by Spencer W. Thomas.
  387.    Note that tran has to have at least 16 entries which will be set
  388.    to:
  389.       Sx, Sy, Sz, Shearxy, Shearxz, Shearyz, Rx, Ry, Rz, Tx, Ty, Tz,
  390.       P(x, y, z, w)
  391. */
  392. int
  393. lib_tx_unwind(tx_mat, tran)
  394.    MATRIX tx_mat;
  395.    double *tran;
  396. {
  397.    int i, j;
  398.    MATRIX locmat, pmat, invpmat;
  399.    COORD4 prhs, psol;
  400.    COORD3 row[3];
  401.  
  402.    lib_copy_matrix(locmat, tx_mat);
  403.  
  404.    /* Divide through by the homogenous value (normalize) */
  405.    if (locmat[3][3] != 0.0)
  406.       for (i=0;i<4;i++)
  407.      for (j=0;j<4;j++)
  408.         locmat[i][j] /= locmat[3][3];
  409.  
  410.    /* pmat is used to solve for perspective, but it also provides
  411.       an easy way to test for singularity of the upper 3x3 component */
  412.    lib_copy_matrix(pmat, locmat);
  413.    for (i=0;i<3;i++)
  414.       pmat[i][3] = 0.0;
  415.    pmat[3][3] = 1.0;
  416.  
  417.    if (lib_matrix_det4x4(pmat) == 0.0)
  418.       return 0;
  419.  
  420.    /* First, isolate perspective */
  421.    if (locmat[0][3] != 0.0 ||
  422.        locmat[1][3] != 0.0 ||
  423.        locmat[2][3] != 0.0) {
  424.       /* prhs is the right hand side of the equation */
  425.       prhs[X] = locmat[0][3];
  426.       prhs[Y] = locmat[1][3];
  427.       prhs[Z] = locmat[2][3];
  428.       prhs[W] = locmat[3][3];
  429.  
  430.       /* Solve the equation by inverting pmat and multiplying
  431.      prhs by the inverse.  */
  432.       lib_invert_matrix(invpmat, pmat);
  433.       lib_transform_coord(psol, prhs, invpmat);
  434.  
  435.       /* Save the perspective information */
  436.       tran[U_PERSPX] = psol[X];
  437.       tran[U_PERSPY] = psol[Y];
  438.       tran[U_PERSPZ] = psol[Z];
  439.       tran[U_PERSPW] = psol[W];
  440.  
  441.       /* Clear the perspective partition */
  442.       locmat[0][3] = 0.0;
  443.       locmat[1][3] = 0.0;
  444.       locmat[2][3] = 0.0;
  445.       locmat[3][3] = 1.0;
  446.       }
  447.    else {
  448.       /* No perspective */
  449.       tran[U_PERSPX] = 0.0;
  450.       tran[U_PERSPY] = 0.0;
  451.       tran[U_PERSPZ] = 0.0;
  452.       tran[U_PERSPW] = 0.0;
  453.       }
  454.  
  455.    /* Pull out the translation */
  456.    for (i=0;i<3;i++) {
  457.       tran[U_TRANSX+i] = locmat[3][i];
  458.       locmat[3][i] = 0.0;
  459.       }
  460.  
  461.    /* Figure out scale and shear */
  462.    for (i=0;i<3;i++) {
  463.       row[i][X] = locmat[i][0];
  464.       row[i][Y] = locmat[i][1];
  465.       row[i][Z] = locmat[i][2];
  466.       }
  467.    /* Compute X scale factor and normalize the first row */
  468.    tran[U_SCALEX] = lib_normalize_vector(row[0]);
  469.  
  470.    /* Compute XY shear factor and make 2nd row orthogonal to 1st */
  471.    tran[U_SHEARXY] = DOT_PRODUCT(row[0], row[1]);
  472.    COMB_COORD(row[1], row[1], row[0], 1.0, -tran[U_SHEARXY])
  473.  
  474.    /* Compute Y scale and normalize 2nd row */
  475.    tran[U_SCALEY] = lib_normalize_vector(row[1]);
  476.    tran[U_SHEARXY] /= tran[U_SCALEY];
  477.  
  478.    /* Compute XZ and YZ shears, orthogonalize 3rd row */
  479.    tran[U_SHEARXZ] = DOT_PRODUCT(row[0], row[2]);
  480.    COMB_COORD(row[2], row[2], row[0], 1.0, -tran[U_SHEARXZ])
  481.    tran[U_SHEARYZ] = DOT_PRODUCT(row[1], row[2]);
  482.    COMB_COORD(row[2], row[2], row[1], 1.0, -tran[U_SHEARYZ])
  483.  
  484.    /* Get Z scale and normalize 3rd row */
  485.    tran[U_SCALEZ] = lib_normalize_vector(row[2]);
  486.    tran[U_SHEARXZ] /= tran[U_SCALEZ];
  487.    tran[U_SHEARYZ] /= tran[U_SCALEZ];
  488.  
  489.    /* At this point, the matrix (in rows[]) is orthonormal.
  490.       Check for a coordinate system flip.  If the determinant
  491.       is -1, then negate thematrix and the scaling factors */
  492.    CROSS(prhs, row[1], row[2])
  493.    if (DOT_PRODUCT(row[0], prhs) < 0.0)
  494.       for (i=0;i<3;i++) {
  495.      tran[U_SCALEX+i] *= -1.0;
  496.      row[i][X] *= -1.0;
  497.      row[i][Y] *= -1.0;
  498.      row[i][Z] *= -1.0;
  499.      }
  500.  
  501.    /* Get the rotations out */
  502.    tran[U_ROTATEY] = asin(-row[0][Z]);
  503.    if (cos(tran[U_ROTATEY]) != 0) {
  504.       if (fabs(row[1][Z]) < EPSILON && fabs(row[2][Z]) < EPSILON)
  505.      tran[U_ROTATEX] = 0.0;
  506.       else
  507.      tran[U_ROTATEX] = atan2(row[1][Z], row[2][Z]);
  508.       if (fabs(row[0][Y]) < EPSILON && fabs(row[0][X]) < EPSILON)
  509.      tran[U_ROTATEZ] = 0.0;
  510.       else
  511.      tran[U_ROTATEZ] = atan2(row[0][Y], row[0][X]);
  512.       }
  513.    else {
  514.       tran[U_ROTATEX] = atan2(row[1][X], row[1][Y]);
  515.       tran[U_ROTATEZ] = 0.0;
  516.       }
  517.  
  518.    return 1;
  519. }
  520.